%%%%%
% CODE DESCRIPTION :
%     load3ds reads in the spectroscopic file generated by the Nanonis
%     software. It parses the header and reads in the Binary data as well.
%     The return value contains an array of structure of different
%     channels(data) and the detail experimental parameters(exp_par) of
%     each coordinates.
%
% CODE HISTORY :
%     06/10/2013 HENRY YU CREATED
%     18/01/2014 Guan add MLS mode
%     27/12/2014 Chih-Chuan Su modify MLS mode for Nanonis of B210
%     19/03/2016 Chih-Chuan Su modify file reading for unfinished data
%     29/03/2016 Chih-Chuan Su read bias of topo if exist
%     19/09/2017 Chih-Chuan Su rewrite the part of reading data to speed up
%%%%%
function [data, exp_par] = load_3ds(filename, pathname, calibration)
header.filename = filename;
header.pathname = pathname;

if exist([pathname filename], 'file')
    fid = fopen([pathname filename], 'r', 'ieee-be');    % open with big-endian
else
    fprintf('File does not exist.\n');
    return;
end

% read header data
    %{
        The header consists of key-value pairs, separated by an equal sign,
        e.g. Grid dim="64 x 64". If the value contains spaces it is enclosed by
        double quotes (").
    %}
    while 1
        s = strtrim(fgetl(fid));
        if strcmpi(s,':HEADER_END:')
            break
        end

        s1 = strsplit(s,'=');
        s_key = strrep(lower(s1{1}), ' ', '_');
        s_val = strjoin(s1(2:end),'=');
        s_val = strrep(s_val, '"', '');
        switch s_key
            case 'grid_dim'
                s_val = strsplit(s_val, 'x');
                header.icols = str2num(s_val{1});
                header.irows = str2num(s_val{2});

            case 'grid_settings'
                s_grid_settings = sscanf(s_val, '%f;%f;%f;%f;%f');
                header.xcenter = s_grid_settings(1); header.ycenter = s_grid_settings(2);
                header.xrange = s_grid_settings(3); header.yrange = s_grid_settings(4);
                header.angle=s_grid_settings(5);
    %           header.x_distmin =  -W/2; header.y_distmin =  -H/2;
    %           header.x_distmax =   W/2; header.y_distmax =   H/2;
                header.xyunit = 'm';

            case 'filetype'
                header.filetype = s_val;

                % get MLS parameters
            case 'segment_start_(v),_segment_end_(v),_settling_(s),_integration_(s),_steps_(xn),_lockin,_init._settling_(s)'
                mls_sets = strsplit(s_val, ';');
                mls_num=size(mls_sets,2);
                for i=1:mls_num
                    mls_scanval = sscanf(mls_sets{:,i}, '%f,%f,%f,%f,%f');
                    header.mls_para.starV(i) = mls_scanval(1);
                    header.mls_para.endV(i) = mls_scanval(2);
                    header.mls_para.settlS(i) =  mls_scanval(3);
                    header.mls_para.intS(i) =  mls_scanval(4);
                    header.mls_para.steps(i) =  mls_scanval(5);
                end

                % fixed parameters, experiment parameters, channels:
            case {'fixed_parameters', 'experiment_parameters', 'channels'}
                s_vals = strsplit(s_val, ';');
                header.(s_key) = s_vals;

                % number of parameters
            case '#_parameters_(4_byte)'
                header.num_parameters = str2num(s_val);

                % experiment size
            case 'experiment_size_(bytes)'
                header.experiment_size = str2num(s_val);

                % spectroscopy points
            case 'points'
                header.ilayers = str2num(s_val);
                %         header.points = str2num(s_val);

                % delay before measuring
            case 'delay_before_measuring_(s)'
                header.delay_before_meas = str2num(s_val);

                %comment
            case 'comment'
                header.comment=s_val;

                % other parameters -> treat as strings
            otherwise
                s_key = regexprep(s_key, '[^a-z0-9_]', '_');
                if length(s_key) > 30
                    s_key=s_key(1:30);
                end
                header.(s_key) = s_val;
        end
    end

%read the raw data
    dim = [header.irows, header.icols];
    num_channels = size(header.channels,2);
    num_data_set = header.num_parameters+(header.ilayers*num_channels);
    num_data_points =  dim(1)*dim(2);
    num_data_all = num_data_points*num_data_set;

    [input_1d, count] = fread(fid,'float');
    if count < num_data_all
        input_1d(count+1:num_data_all) = nan;
    end

%reshpae the input
    input_2d=reshape(input_1d,[num_data_set,num_data_points])';

    input_2d_par=input_2d(:,1:header.num_parameters);

    input_2d_exp=input_2d(:,header.num_parameters+1:num_data_set);
    input_2d_exp = reshape(input_2d_exp, dim(1), dim(2), header.ilayers, num_channels);
% collect the important parameters
    % file name
    name = strrep(filename(1:end-4), ' ', '');
    
    %creat 'exp_par' file
    exp_par.name = LegalizeName(name);
    exp_par.map = read_exp_par(input_2d_par, header);
    exp_par.map = flipud(exp_par.map);

    % energy list
    if isfield(header,'filetype')
        if strcmp(header.filetype, 'MLS')
            % MLS energy list
            e=[];
            for i=1:mls_num
                e = [e linspace(header.mls_para.starV(i),header.mls_para.endV(i),header.mls_para.steps(i))];
            end
            e=unique(e, 'stable');
        end
    end
    if ~isfield(header,'filetype') || strcmp(header.filetype, 'Linear')
        sweep_start = exp_par.map.sweep_start(1,1);  %header.info.bias_spectroscopy_sweep_start__v_
        sweep_end = exp_par.map.sweep_end(1,1);  %header.info.bias_spectroscopy_sweep_end__v_
        e = linspace(sweep_start, sweep_end, header.ilayers);
    end

    % the x, y coordinate matrix
    [rx,ry]=meshgrid(linspace(0,header.xrange,dim(2)+1), linspace(header.yrange,0,dim(1)+1));
    rx=(rx(1:dim(1),1:dim(2)) + rx(1:dim(1),2:dim(2)+1))/2;
    ry=(ry(1:dim(1),1:dim(2)) + ry(2:dim(1)+1,1:dim(2)))/2;
    [Rx_tmp, Ry_tmp] = meshgrid(linspace(-header.xrange/2 ,header.xrange/2, dim(2)+1), linspace(header.yrange/2, -header.yrange/2, dim(1)+1));
    Rx_tmp=(Rx_tmp(1:dim(1),1:dim(2)) + Rx_tmp(1:dim(1),2:dim(2)+1))/2;
    Ry_tmp=(Ry_tmp(1:dim(1),1:dim(2)) + Ry_tmp(2:dim(1)+1,1:dim(2)))/2;
    Rx=Rx_tmp*cosd(header.angle)+Ry_tmp*sind(header.angle)+header.xcenter;
    Ry=-1*Rx_tmp*sind(header.angle)+Ry_tmp*cosd(header.angle)+header.ycenter;

%create maps files for different channel
data(num_channels).map = zeros(dim(1), dim(2), header.ilayers);
for i = 1:num_channels
    data(i).type = '3D_map';
    if ~isnan(calibration)
        data(i).map = flip(input_2d_exp(:,:,:,i),1)*calibration*1e9;
    else
        data(i).map = flip(input_2d_exp(:,:,:,i),1);
    end
    data(i).map = permute(data(i).map , [2,1,3]);
    data(i).map = flip( flip( data(i).map, 2) , 1);
    channel_name = strrep(header.channels{i}, ' ', '');
    unit_name = channel_name(end-2:end);
    if strcmpi(unit_name, '(A)')
         if ~isnan(calibration)
            data(i).map_unit = 'nS';
         else
             data(i).map_unit = 'A';
         end
    elseif strcmpi(unit_name, '(m)')
        data(i).map_unit = 'm';
    else
        data(i).map_unit = 'unit';
    end

    data(i).channel = channel_name(1:end-3);
    data(i).name = LegalizeName(name);
    data(i).e = e;
    data(i).rx = rx;
    data(i).ry = ry;
    data(i).Rx = Rx;
    data(i).Ry = Ry;
    data(i).info = header;
    data(i).sweep_unit = header.sweep_signal(strfind(header.sweep_signal,'(')+1 : strfind(header.sweep_signal,')')-1);
    %DataProInfo:
    data(i).DataProInfo{1}=0;
    data(i).DataProInfo{2}{1,1}='raw';data(i).DataProInfo{2}{1,3}='';
    data(i).DataProInfo{2}{1,2}=[data(i).name,'_',ChannelAbbr(data(i).channel)]; 
    data(i).DataProInfo{3}=zeros(header.ilayers,2);
    for j=header.ilayers:-1:1
        tmp=data(i).map(:,:,j);
        m = nanmean(tmp(:));
        s = nanstd(tmp(:));        
        data(i).DataProInfo{3}(j,:)=[m,s];
        data(i).DataProInfo{5}(j,:)=[-4,4];
    end
    data(i).DataProInfo{4}=2;
    
    data(i).DataProInfo{6}=[header.xcenter header.ycenter,header.xrange,header.yrange,...
        -header.angle,header.icols,header.irows];
    
    %replace nan with 0
    data(i).map(isnan(data(i).map))=0;
end

% add an additional section for the Topograph parameters
topo_map = flipud(exp_par.map.z_m);
data(num_channels+1) = data(num_channels);
data(num_channels+1).map = topo_map;
data(num_channels+1).channel = 'Z';
data(num_channels+1).map_unit = 'm';
if isfield(header,'bias_bias__v_')
    data(num_channels+1).e = str2double(header.bias_bias__v_);
else
    data(num_channels+1).e = 0;
end
%DataProInfo:
data(num_channels+1).DataProInfo{1}=0;
data(num_channels+1).DataProInfo{2}{1,1}='raw';data(num_channels+1).DataProInfo{2}{1,3}='';
data(num_channels+1).DataProInfo{2}{1,2}=[data(num_channels+1).name,'_',ChannelAbbr(data(num_channels+1).channel)];
data(num_channels+1).DataProInfo{3}=...
    [nanmean(data(num_channels+1).map(:)),nanstd(data(num_channels+1).map(:))];
data(num_channels+1).DataProInfo{4}=2;
data(num_channels+1).DataProInfo{5}(1,:)=[-4,4];
data(num_channels+1).DataProInfo{6}=[header.xcenter header.ycenter,header.xrange,header.yrange,...
    -header.angle,header.icols,header.irows];

%replace nan with 0
data(num_channels+1).map(isnan(data(num_channels+1).map))=0;

fclose(fid);



%{
    read_exp_par() is a function that reads the experiment parameters of the
	current experiment, i.e., the current grid point. It returns a struct of
    all the parameters of this point.
%}
    function ExpPar = read_exp_par(ParVal, Header)
        ParVal = reshape(ParVal, header.icols, header.irows, header.num_parameters);
        ParVal = permute(ParVal,[2,1,3]);
        par_names = [Header.fixed_parameters Header.experiment_parameters];
        for ii = 1:length(par_names)
            par_name = lower(LegalizeName(par_names{ii}));
            ExpPar.(par_name) = ParVal(:,:,ii);
        end
    end

end